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The Floquet spectra of a class of driven SU(2) systems have been shown to display butterfly 
patterns with multifractal properties. The implication of such critical spectral behavior for the 
Floquet eigenstate statistics is studied in this work. Following the methodologies for understanding 
the fractal behavior of energy eigenstates of time-independent systems on the Anderson transition 
point, we analyze the distribution profile, the mean value, and the variance of the logarithm of 
the inverse participation ratio of the Floquet eigenstates associated with multifractal Floquet spec- 
tra. The results show that the Floquet eigenstates also display fractal behavior, but with features 
markedly different from those in time-independent Anderson-transition models. This motivated 
us to propose a new type of random unitary matrix ensemble, called "power-law random banded 
unitary matrix" ensemble, to illuminate the Floquet eigenstate statistics of critical driven systems. 
The results based on the proposed random matrix model are consistent with those obtained from 
our dynamical examples with or without time-reversal symmetry. 

PACS numbers: 05.45. Df, 05.45. Mt, 71.30.+h, 74.40. Kb, 05.45.-a 



I. INTRODUCTION 

The critical behavior of time- independent systems, es- 
pecially in terms of the spectral statistics and the eigen- 
state statistics, has attracted great attention. On the 
spectrum side, Hofstadter's butterfly spectrum of the 
Harper model has been a paradigm for critical spec- 
tral statistics, representing a multifractal spectrum [l|, Q 
of a system exactly on the metal-insulator transition 
point. On the eigenstate side, mainly through studies in 
time-independent models, such as the power-law random 
banded matrix (PRBM) model [|[ and the standard An- 
derson tight-binding model (TBM) it has been well- 
established that for a system on a metal-insulator transi- 
tion point or the Anderson transition point [5] , its eigen- 
states show clear fractal features. This background of 
understanding the critical behavior of time-independent 
systems motivated our interest in the critical behavior of 
periodically driven systems. Below we first introduce re- 
cent related studies of critical Floquet spectra, and then 
briefly describe the motivation and the results of this 
work. 

It is well known that the Floquet (quasi-energy) spec- 
trum of a delta-kicked version of the Harper model also 
displays Hofstadter's butterfly patterns jaQ- Interest- 
ingly, though the kicked Harper model (KHM) can be 
classically chaotic, its spectrum, due to its fractal na- 
ture, does not follow the Bohigas-Giannoni-Schmit con- 
jecture [|[ at all. This makes the KHM not only a fruitful 
model for gaining new insights into the issue of quantum- 
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classical correspondence in classically chaotic systems, 
but also an intriguing model to study critical spectral 
statistics. Indeed, for quite a long time, studies of frac- 
tal Floquet spectra were largely restricted to the KHM 
and its variants . In a proposal to experimentally real- 
ize Hofstadter's butterfly Floquet spectrum in cold-atom 
laboratories, Wang and Gong [HI [H| recently demon- 
strated that Hofstadter's butterfly Floquet spectrum can 
be synthesized by use of a double-kicked cold-atom ro- 
tor system [l2j under a quantum resonance condition. 
Lawton et al. [l3[ then showed that the butterfly Flo- 
quet spectrum of the cold-atom system studied in Refs. 
[ru [UJ is equivalent to that of the standard KHM if and 
only if one system parameter takes irrational values. In 
addition to motivating a cold-atom realization of criti- 
cal Floquet spectra of periodically driven systems, Refs. 
[Tfl [TlT | seem to have offered a general strategy for syn- 
thesizing critical Floquet spectra in driven systems. 

Using an approach extended from Refs. [lol . fllj |. re- 
cently Wang and Gong [l4[ showed that the Floquet spec- 
tra of a class of driven SU(2) systems also display but- 
terfly patterns and multifractal properties that are char- 
acteristics of highly critical spectra. This establishes a 
completely different class of critical driven systems with- 
out a connection with the KHM context. Interestingly, 
the driven SU(2) model in Ref. T4] can be understood 
as a simple extension of the well-known kicked top model 
(KTM) [l5[ in the quantum chaos literature. Because 
the KTM has just been experimentally realized in a cold 
133 Cs system [Ig, it can be expected that a critical driven 
SU(2) system may also be experimentally realized using 
the collective spin of a 133 Cs atomic ensemble. An alter- 
native experimental realization may be based on a driven 
two-mode Bose- Einstein condensate [13, 53, Gil, which 



2 



represents a strongly self-interacting driven system. 

Given the above-mentioned class of driven quantum 
systems with critical Floquet spectra, it becomes nec- 
essary to study the behavior of the associated Flo- 
quet eigenstates. Theoretically speaking, because driven 
SU(2) systems always have a finite number of Floquet 
eigenstates, the eigenstate analysis becomes much easier 
than in the KHM, with the latter necessarily involving an 
infinite number of eigenstates for a fractal Floquet spec- 
trum. A careful investigation of the Floquet eigenstates 
over the entire spectrum will help to better understand 
the critical behavior in time-dependent systems in gen- 
eral. Experimentally speaking, information about the 
eigenstate statistics may be more directly accessible to 
measurements than a fractal spectrum. 

To analyze the critical behavior of the Floquet eigen- 
states in driven SU(2) systems, we adopt the same ap- 
proach as in previous studies of time-independent sys- 
tems. That is, we shall numerically examine the fluctua- 
tions of the eigenstates pjj]. The eigenstate fluctuations 
can be characterized by a set of inverse participation ra- 
tios (IPR): 
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where A is the index of the eigenstates, \cj>\) represents 
one eigenstate under investigation, and {\n}} are the 
basis states. For convenience we focus on the IPR P2 
(i.e., q = 2). By analogy to critical eigenstate behavior 
in time-independent systems, we expect that P2 scales 
anomalously with the Hilbert space dimension N as 
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where is a fractal dimension of a particular eigen- 
state 1 4>\). But is there also a unique fractal dimension 
D 2 for the average behavior of all the Floquet eigenstates, 
for example, via the slope of the averaged ln(P 2 ), denoted 
(ln(P2)), versus ln(iV)? To that end, we shall examine 
if, as the system gets closer to the thermodynamic limit 
(N — » +00), the distribution of ln(P 2 ) shows signs of a 
scale-invariant form [2(|. In other words, whether the 
distribution function of ln(P2), denoted II[ln(P2)], only 
shifts as N varies. 

Certainly, the system under our study has only a finite 
size N. In time-independent Anderson-transition studies 
using the PRBM or the TBM, it was conjectured that 
the variance of ln(P2), denoted ex 2 (TV), scales with N as 
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with a 1 (00), A, and 7 being three adjustable parame- 
ters [2lJ. For a d-dimcnsional system on the Anderson 
transition point, it was shown that 7 is related to D% by 
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where /3 equals 1 or 2 depending upon whether or not 
the system has time- reversal symmetry 22]. As one 
main task of this work, we shall examine if these re- 
sults for time-independent systems still hold for critical 
Floquet eigenstates. Furthermore, we hope to see how 
the criticality of the eigenstates of unitary operators dif- 
fers from the criticality of the eigenstates of self-adjoint 
operators. Results along this direction will also be rel- 
evant to recent investigations on the "unitary Anderson 
model" [23|, the Thue-Morse sequence generating multi- 
fractal eigenstates of the quantum baker's map 24[, the 
one-parameter model of quantum maps showing multi- 
fractal eigenstates 25| . as well as recent experimental 
and theoretical studies of Anderson transition in kicked- 
rotor systems [26|, (2?J ■ 

We now briefly summarize the main findings of this 
work. For the driven SU(2) systems studied here, we 
consider two different parameter regimes: in one regime 
the Floquet spectra display clear butterfly patterns, and 
in the other regime, the butterfly patterns of the Floquet 
spectra have dissolved due to increased strength of the 
driving fields. For both regimes, we find that II[ln(P 2 )] 
is not as smooth as observed in the TBM or PRBM, 
indicating some non-universal features in dynamical sys- 
tems. The n[ln(P 2 )] for cases with dissolved butterfly 
patterns is however smoother. For either regime, it is 
found that the ensemble average (ln(p2)) does scale lin- 
early with ln(iV), with the slope of the (ln(P2)) vs ln(iV) 
curve clearly defining the fractal dimension D 2 for all the 
eigenstates. We also find it possible to fit the variance of 
ln(P 2 ) by Eq. ([3]), but with the exponent 7 given by 
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instead (with d = 1); i.e., a factor of two is missing 
from the denominator as compared with Eq. (4) for 
time-independent critical systems. To further under- 
stand this difference, we propose a random matrix model, 
which we call "power-law random banded unitary ma- 
trix" (PRBUM) model. By tuning the parameters of the 
PRBUM, the D 2 value associated with the PRBUM can 
be varied. More interestingly, we observe that the vari- 
ance a 2 (N) of the PRBUM also follows Eq. $S§, with the 
exponent 7 again given by Eq. §5§. This suggests that 
our findings about the Floquet eigenstate statistics based 
on driven SU(2) systems do reflect some general aspects 
of critical Floquet eigenstates. 

This paper is organized as follows. In Sec. [IIJ we 
present detailed results of the eigenstates statistics in 
our driven SU(2) models, with or without time-reversal 
symmetry. In Sec. IIII1 we introduce the PRBUM to rep- 
resent a class of critical Floquet operators, discuss the 
statistics of the eigenstates of PRBUM, and then com- 
pare the associated results with those found in actual 
dynamical systems. In Sec. IIV[ we study the eig enstate 
statistics of the standard kicked top model 15|, which 



represents a classically chaotic, but non-critical, driven 
system. Section IVl concludes this work. 
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II. FRACTAL STATISTICS OF THE FLOQUET 
EIGENSTATES IN DRIVEN SU(2) MODELS 

The focus in Ref. 14] is on the fractal spectral statis- 
tics. Here, using the same model we study the statistics 
of the Floquet eigenstates. The first Floquet operator 
under study is given by 



exp(~iaJ x ) exp y^^-^JJ exp(— ia-J x ), 

(6) 

where J Xl J y , J z are angular momentum operators satis- 
fying the SU(2) algebra and J is the conserved total angu- 
lar momentum quantum number that defines a (2 J + 1)- 
dimensional Hilbert space. Readers can refer to Ref. [14| 
for detailed descriptions and motivations of this model. 
This model is also called as a "double-kicked top model" 
(DKTM) in Ref. [lj]. 

Eigenstates of the J z operator are denoted as |m), with 
Jz\ m ) = m\m). States will be chosen as our repre- 

sentation for eigenstate analysis. To analyze the Floquet 
eigenstates, it is necessary to express the Floquet oper- 
ator in symmetric basis states, a procedure that block- 
diagonalizes the Floquet matrix. On the one hand, this 
will simplify our analysis; on the other hand, this is nec- 
essary for the sake of comparison between an actual dy- 
namical system and the PRBUM model proposed below. 

The DKTM Floquet operator F in Eq. JBJ has a parity 
unitary symmetry R*<FR = F where R = exp(-iirJ x ). 
This symmetry can be used to block diagonalize the 
F matrix into two disconnected sub-matrices associated 
with either odd-parity or even-parity subspaces. With- 
out loss of generality we only present below results for the 
J-dimensional odd-parity subspace. Besides the parity 
symmetry, F also has a time-reversal anti-unitary sym- 
metry TFT = F+, with 
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where K is the complex conjugation operator. To explore 
the implication of this time-reversal symmetry for the 
eigenstate statistics, we shall also consider a variant of 
F, i.e., 

/ rjJ 2 \ ( r)J 2 \ 

F' = exp f i-^j J exp(—iaJ x ) exp f -i-^j J exp(-ia J y ). 

' ' ' "' ' (8) 

Evidently, F' differs from F only in the last factor, i.e., 
exp(— iaJ x ) in F is replaced by exp(— ia J y ). Because of 
this difference, we call F in Eq. (j6]) the J x — J x model 
and call F' the J x — J y model. It is easy to check that 
the latter does not have the parity symmetry or the time- 
reversal symmetry. For the J x — J y model, which cannot 
be reduced to any block-diagonal form, we examine the 
eigenstates of the full Floquet matrix. 

For both cases we define a dimensionless system pa- 
rameter, 
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This choice of being tt times the golden mean is to 
ensure that the resulting Floquet eigenstate statistics is 
indeed representative of driven systems with fractal Flo- 
quet spectra. As detailed below, we consider two dif- 
ferent regimes for the product aJ. In the first regime 
defined by 0.95 < aJ < 1.05, the Floquet spectra show 
clear butterfly patterns; in the second regime defined by 
9.95 < aJ < 10.05, the butterfly spectra have dissolved, 
with fractal dimensions of the spectra increased [l4[ ■ 



A. J x — Jx model 

This is a time-reversal symmetric system. Because 
Dyson's circular ensemble of random unitary matrices 
|15j with time-reversal symmetry is called "circular- 
orthogonal-ensemble" (COE), we regard the J x — J x 
model as an example of critical COE statistics. 



1. 0.95 <aJ < 1.05 

Figure [lj a) shows the distributions of the logarithm of 
the IPR P2, denoted II[ln(P2)], for different J. It is seen 
that the distribution function n[ln(P2)] is not as smooth 
as that observed in early Anderson-transition studies [20I — 
I22I ]. Nevertheless, it is clear that as J increases, the left 
tail of n[ln(P2)] systematically shifts to the left direc- 
tion associated with more negative ln(P 2 ). The profile 
of II[ln(P2)], though somewhat changes as J increases, 
does maintain its main features as J increases. Due to 
these features that are similar to early findings for the 
critical eigenstates in time-independent systems, it can 
be expected that the average of ln(P 2 ) will show a scal- 
ing behavior with ln(J). As shown in Fig. [ljb), this is 
indeed the case. Therein, (ln(P2)), obtained by averaging 
ln(P 2 ) over all eigenstates (in the odd-parity subspace), 
displays an excellent linear behavior with ln(J). From 
the slope of the fitting line in Fig. 1(b), we are able to 
obtain the fractal dimension D 2 — 0.274. 

The distribution profile n[ln(P2)] in Fig. 1(a) is seen 
to display rich features, with significant fluctuations and 
multiple notable peaks. Qualitatively, this reflects that 
our system is an actual dynamical system and hence 
the underlying rich dynamics will manifest itself through 
some non- universal statistical features. Related to this 
observation we also note that in our calculations, all 
the Floquet eigenstates are treated equally and all of 
them are used for averaging. This is in contrast to the 
common procedure in analyzing time-independent crit- 
ical systems, where only those energy eigenstates in a 
certain small energy window around zero eigenvalue are 
included to examine the distribution of ln(P 2 ) [2(]|-[22]. 
The justification for including all Floquet states in our 
analysis is as follows: the quasi-energy spectra lie on a 
unit circle and hence all states with different eigenphases 
on the unit circle should be treated on equal footing. To 
double check this understanding, we have also taken win- 
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dows of different widths centered around zero value of the 
eigenphase and then calculate the distribution of ln(P 2 ). 
No improvement in the smoothness of II(lnP 2 ) is found. 
Rather, we obtained similar distribution of ln(P 2 ) with 
clear fluctuations. It is also tempting to connect the non- 
universal features of II(lnP2) with the phase space struc- 
tures of the underlying classical limit. However, such 
a perspective, which calls for a good understanding of 
quantum-classical correspondence in critical systems, is 
unlikely to succeed because the classical limit of our dy- 
namical model is completely chaotic [14J. 

In Fig. Die), we plot ln[cr 2 (oo) - a 2 (J)} vs In (J) (filled 
circles), where cr 2 (J) is the variance of ln(P 2 ) and a 2 (oo) 
is a fitting parameter, whose value is found by fitting our 
data points with the empirical formula given in Eq. ([3]). 
As seen in Fig.[ljc), the fitting is reasonably good, yield- 
ing that [cr 2 (co) — (T 2 (J)] scales as J -7 , with 7 = D 2 
[P/ 2 is already determined by the fitting in Fig. H£b)], 
cr 2 (oo) ~ 0.68, and A ~ 1.40. Despite obvious fluctua- 
tions around the fitting curve, the result in Fig. QJc) sug- 
gests that the tool borrowed from traditional Anderson- 
transition studies for time-independent systems can be 
still useful here. Furthermore (probably more interest- 
ingly), the fitting in Fig. QJc) also unexpectedly reveals 
a big difference from what can be expected from Eq. 
(|4| with d = 1 and j3 = V. here 7 = D 2 instead of 
D2/2. Therefore, an intriguing difference between time- 
independent critical systems and periodically driven crit- 
ical systems is observed here. 



2. 9.95 <aJ < 10.05 

As mentioned above, for this parameter regime the 
butterfly patterns in the Floquet spectra have dissolved 
almost completely. We present the associated eigenstate 
statistics in Fig. [2J In Fig. [2fa), we show the distri- 
bution profile of ln(P2) for different J. In contrast to 
the previous case shown in Fig. 1(a), n[ln(P2)] is now 
much smoother (essentially only one peak is left). From 
the same panel, we also see a systematic left-shift of the 
distribution function as J increases. This systematic left- 
shift leads to an evident linear behavior of the average 
value of ln(P 2 ) as a function of ln( J), as shown in Fig. 
[2jb). The slope of the fitting line in Fig. [2{b) gives 
the fractal dimension D2 — 0.256. Comparing this re- 
sult with that in Fig. 1(b), one sees that though the 
fractal dimension of the Floquet spectra increases due to 
increasing a J [l4| , the fractal dimension of the associated 
eigenstates may decrease. 

In Fig. [2][c), we examine the variance of ln(P 2 ) as a 
function of ln(J) (again, for the odd-parity subspace). 
Same as in Fig. 1(c), we fit our results with the empirical 
formula given in Eq. ©. The fitting in Fig. OJc) is better 
than that in Fig. 1(c), consistent with the fact that the 
distribution of ln(P 2 ) is quite smooth here. The fitting 
in Fig. dc) gives cr 2 (oo) ~ 0.77, A ~ 1.59, and 7 = D 2 , 
where the value of Z? 2 is found in Fig. HJb) . The finding 
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FIG. 1: (Color online) (a) Distribution of ln(P2) for the 
J x - J x model, with J = 200, 400, 800, 1600 and 3200, in 
the representation of odd-parity basis states denned in the 
text. The size of the Floquet matrix ensemble is important 
for numerical simulation. In order to construct the necessary 
ensemble, we consider a range of a, i.e., 0.95 < aj < 1.05, 
yielding respectively 4000, 2000, 1000, 500 and 250 matrices 
for the different values J. (b) (ln(P2)}, the mean value of 
ln(P2) averaged over all Floquet eigenstates, as a function of 
ln( J). The slope of the fitting line gives D2 — 0.274. (c) Loga- 
rithm of [a 2 (00) — a 2 (J)] as a function of ln(J), where J is the 
dimension of the odd-parity Hilbert subspace. Filled circles 
are our numerical results for the J x — J x model, and the solid 
line is the fitting of the numerical results using the empirical 
formula given in Eq. © with a 2 {oo) = 0.68, A = 1.40 and 
7 = Di. The plotted variables here and in all other figures 
are dimensionless. 



that 7 is not equal to P/ 2 /2 but D 2 again strengthens our 
early observation from Fig. 1. 



B. J x 



To verify if our findings above are general, we now 



turn to the J x — J y model [Eq. 



Due to the lack of 



time-reversal symmetry here, this case can be regarded as 
an example of critical "circular-unitary-ensemble" (CUE) 
statistics. All the eigenstates of the Floquet operator F' 
will be considered. 



1. 0.95 < aJ < 1.05 

For this regime where the butterfly patterns of the 
Floquet spectra can be clearly seen, Fig. [3fa) displays 
the distribution of ln(P 2 ) for different Hilbert space di- 
mension N = 2 J + 1. Analogous to the previous case 
with time-reversal symmetry, II[ln(P 2 )] displays interest- 
ing fluctuations. As N increases, II[ln(P 2 )] undergoes 
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In (J) In (J) 

FIG. 2: (Color online) (a) Distribution of ln(p2) for the J x — 
J x model with 9.95 < aj < 10.05 is presented for the odd- 
parity subspace. Other parameters are the same as in Fig. [1] 
(b) (ln(P2)) is plotted as a function of ln(J). The slope of 
the fitting line gives D 2 ~ 0.239. (c) Logarithm of [a 2 (oo) - 
a 2 (J)] as a function of ln( J), where J is the dimension of the 
odd-parity Hilbert subspace. Filled circles are our numerical 
results for the J x — Jx model, and the solid line is the fitting 
of the numerical results using the empirical formula given in 
Eq. ©, with cr 2 (oo) ~ 0.77, A = 1.59 and 7 = D 2 . 



changes in its profile, shifts its left tail, but also main- 
tains many features. In Fig. [3(b) we obtain again a nice 
linear scaling behavior of (ln(p2)) with ln(iV). From the 
slope of the linear scaling, we obtain the fractal dimen- 
sion D 2 — 0.259. This D 2 value is different from that for 
the J x — J x model with the same values of a J(Note that 
the spectral statistics for the J x — J y model also differs 
from that for the J x — J x model 14]). 

Same as in Fig. QJc), in Fig. |3Jc) we study the 
variance of ln(p2) [now denoted a 2 (N)] as a function 
of ln(iV), using the fitting formula given in Eq. (J3|). 
The fitting, though with clear fluctuations, yields that 
[cr 2 (oo) - a 2 {N)] scales as N~i, with cr 2 (oo) ~ 0.92, 
A ~ 1.04, and 7 = D 2 /2 [D 2 value obtained from Fig. 
|3Jb)]. Remarkably, though Eq. (4) with d = 1 and (3 = 2 
(because of the lack of time-reversal symmetry) predicts 
7 = -D2/4, here we have 7 = D 2 /2 instead. The impor- 
tant common feature shared by the J x — J y model and 
the J x — J x model is thus the missing of a factor of 2 in 
the numerically obtained 7 value as compared with the 
empirical formula for time-independent critical systems. 
This interesting finding also supports the use of Eq. d3j 
as a tool for understanding Floquet eigenstate statistics. 
Our numerical observations here will be further strength- 
ened by a random matrix model. 



2. 9.95 < aJ < 10.05 

Just like the J x — J x model, in this regime the butterfly 
patterns of the Floquet spectra have dissolved. The sta- 
tistical properties of the Floquet eigenstates are shown in 
Fig.|H In Fig.@|a), the distributions of ln(p2) is seen to 
be much smoother than those seen in Fig. 3(a). This is 
somewhat expected from our early findings in the J x — J x 
model. Figure 2(b) shows a linear scaling of (ln(P2)) vs 
ln(iV), with its slope giving D 2 ~ 0.177. In Fig. IUc), we 
study the variance of \n(P 2 ) as a function of ln(iV), as 
compared with the empirical formula given in Eq. ([3j: 
the fitting with the empirical formula is excellent, yield- 
ing cr 2 (oo) ~ 1.11, A ~ 1.17, and 7 = D 2 /2, where the 
value of D 2 is determined in Fig. QJb). Once again, here 
we find 7 = D 2 /2 instead of 7 = D 2 /4: [as suggested by 
Eq. (4) with = 2]. 



III. EIGENSTATE STATISTICS OF PRBUM 

In studies of time-independent critical systems, the 
PRBM model at criticality Q has proved to be fruitful. 
The PRBM is an ensemble of random Hermitian matri- 
ces whose matrix elements {-ffy } are independently dis- 
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FIG. 3: (Color online) (a) Distribution of ln(P2) for the J x — J y 
model, with J = 100(4000), 200(2000), 400(1000), 800(500), 
and 1600(250). The numbers in the brackets are the size of the 
Floquet matrix ensemble. The dimension of the Hilbert space 
is given by N = 2,7 + 1. (b) (ln(P2)), the ensemble mean value 
of ln(P2) averaged over all Floquet eigenstates, as a function 
of ln(JV). The slope of the curve of (ln(P2)) vs ln(iV) gives 
D 2 ~ 0.259. (c) Logarithm of [cr 2 (oo) - a 2 (N)} vs ln(JV). 
Filled circles are numerical results for the J x — J y model, 
and the solid line is the fitting of the numerical results using 
the empirical formula given in Eq. ©, with a 2 (00) = 0.92, 
A = 1.04, and 7 = D 2 /2. 
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tributed Gaussian random numbers with mean (-ffy) = 
and the variance satisfying 



<j 2 (H l3 ) 



*-3\ 
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(10) 



The case g = 1 represents the critical point and < 
b <C oo is a parameter characterizing the ensemble. A 
straightforward interpretation of this model is that it 
describes a one-dimensional sample with random long- 
range hopping, with the hopping amplitude decaying as 
\i — j Motivated by our results above for critical Flo- 
quet states, we aim to propose a class of random unitary 
matrices, whose Floquet eigenstate statistics can show 
some general aspects of critical statistics and can be used 
to shed some light on actual dynamical systems. Our nat- 
ural starting point for generating such random unitary 
matrices are the Hcrmitian PRBM. 



A. Algorithm 

To generate a random unitary matrix from a Hermitian 
matrix in the PRBM ensemble, we employ the algorithm 
by Mezzadri, whose original motivation is to generate 
CUE random matrices 28[ from general Gaussian ran- 
dom matrices. For the sake of completeness, we have 
presented a description of Mezzadri's algorithm in Ap- 
pendix [A] For our purpose, that is, to generate a critical 




FIG. 4: (Color online) (a) Distribution of ln(P 2 ) for the J x - 
J y model with 9.95 < aj < 10.05 is presented for different 
Hilbert space dimension N = 2 J + 1. Other parameters are 
the same as in Fig. [3] (b) (ln(P2)) is plotted as a function 
of ln(iV). The slope of this linear curve gives D2 — 0.177. 
(c) Logarithm of [<r 2 (oo) — <r 2 (iV)] vs ln(iV). Filled circles are 
our numerical results for the J x — J„ model, and the solid 



random unitary matrix, we first set the starting point of 
Mezzadri's algorithm as a PRBM ensemble at the crit- 
ical point (g = 1.0). We then generate an ensemble of 
random unitary matrices (denoted U) of the CUE class. 
Significantly, because of the use of PRBM as the input 
for Mezzadri's algorithm, we find that the variance of 
the matrix elements {Uij} thus obtained also satisfies a 
power-law, i.e., 



J\ 



2<?o* 



(11) 



Here the parameter eto is a common prefactor of the ma- 
trix elements, which can be determined by the unitary 
condition. The parameters go and 60 in Eq. (TTTI) depend 
on the parameters g and b of the PRBM used. As three 
computational examples, panels (d)-(f) of Fig. [5]present 
the dependence of \n[ao / a 2 (Uij) — 1] upon In \i — j\, for 
three ensembles of random unitary matrices we gener- 
ated, with sizes N = 500, 1000, and 2000. If the scaling 
of a 2 (Uij) with is indeed a power law as described 

by Eq. (fTTj) . then one should see a linear dependence of 
ln[ao / (J 2 (Uj) — 1] in ln|i — j\. This is indeed the case 
in Figs. [51(d)- (f). Note that the deviations in Figs. EJd)- 
(f) from the fitting straight lines at very large values of 
In \ i—j\ are due to two trivial reasons. First, for very large 
\i— j\, the value of a 2 (Uij) is vanishingly small and hence 
ln[l / a 2 (Uij) — 1] becomes extremely large, thus yielding 
large fluctuations. Second and more importantly, for a 
fixed matrix size, if \i — j\ is very large, then the avail- 
able number of matrix elements become insufficient for 
good statistics. Indeed, as the matrix size increases from 
N = 500 to N = 2000, it is seen from Figs. E^d)-(f) 
that the validity window of the linear fitting gradually 
extends to larger values of In | i — j \ . 

The random unitary matrices generated in the above 
manner, with their matrix elements satisfying the power- 
law scaling of Eq. (ITT|) . are defined as "power-law ran- 
dom banded unitary matrix" of the CUE type (PRBUM- 
CUE). As detailed in Appendix A, one can then gen- 
erate PRBUM of the COE type (PRBUM-COE) via 
V = UU T . As shown in panels (a)-(c) of Fig. [5j the 
variance of the matrix elements of PRBUM-COE also 
obeys Eq. (fTTj) , with different values of go and b . 

To check whether the PRBUM-COE and PRBUM- 
CUE ensembles show critical statistics, we analyzed their 
eigenstates, especially in terms of the distribution and 
the scaling of ln(P~2)- It is found that as we tune the 
parameter b of the PRBM used in the algorithm, the 
resulting fractal dimensions D 2 can be also tuned con- 
tinuously. For example, the D 2 value of PRBUM can be 
made close to that of our driven SU(2) models. In partic- 
ular, at b = 0.1, we obtain g ~ 0.92 for PRBUM-COE 
and go ^ 0.88 for PRBUM-CUE, yielding D 2 ~ 0.279 
and D2 — 0.251, respectively. These two D2 values are 



line is the fitting of the numerical results using the empirical qu j te close to the D 2 values of the J x - J x and J x - J y 



formula given in Eq. 
7 = D2/2. 



©, with <r 2 (oo) ~ 1.11, A = 1.17 and 



models found in Fig. 1 and Fig. 
these findings in detail. 
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FIG. 5: (Color online) The variance of the matrix elements 
of the random unitary matrices generated by Mazzadri's al- 
gorithm with PRBM as the input. To demonstrate the power 
law scaling, the dependence of ln(ao / a 2 {Uij) — 1) on In \i — j\ 
is plotted, where Uij represents a matrix element at the ith 
row and jth column. Here the PRBM ensemble as the in- 
put is set at the critical point g — 1.0 with the parameter 
b = 0.1. Panels (a)-(c) are for PRBUM-COE, with the di- 
mension N = 500, 1000, and 2000, respectively. The fitting 
function of Eq. |TT]) (solid lines) gives a = 0.165 ± 0.004, 
b = 0.575 ± 0.003, and go = 0.875 ± 0.002. Panels (d)-(f) 
are for PRBUM-CUE, with the dimension N = 500, 1000, 
and 2000, respectively. The fitting function of Eq. (Rip (solid 
lines) yields a = 0.277 ± 0.005, b = 0.355 ± 0.004, and 
go = 0.921 ±0.001. 



PRBUM-COE 




ln(N) 



FIG. 6: (Color online) (a) Distribution of ln(P2) ob- 
tained for PRBUM-COE, with the matrix dimension N = 
200(4000), 400(2000), 800(1000), 1600(500) and 3200(250). 
The numbers in the brackets give the size of the ensemble. 

(b) Same as in Fig. □]» and Fig. 2(b), yielding D 2 ~ 0.279. 

(c) Same as in Fig. [He) and Fig. 2(c), but with J replaced 
by N. The fitting curve gives <t 2 (oo) ~ 0.60, A ~ 1.33, and 
7 = D 2 . 



PRBUM-COE ensemble is around 0.24, which is close to 
the D 2 value previously found in the J x — J x model with 
9.95 < aJ < 10.05. These results clearly support our 
use of PRBUM-COE to illuminate the critical eigenstate 
statistics in the J x — J x model. 



This random unitary matrix ensemble is intended to 
model a critical Floquet operator with time-reversal sym- 
metry. The results for PRBUM-COE generated from 
PRBM with b = 0.1 are shown in Fig. U In Fig. Ua), 
we show the distributions of ln(P2) for different values 
of the matrix dimension N (which is the counterpart of 
J in the J x — J x model) , with all the eigenstates of the 
PRBUM-COE ensemble considered. In contrast to the 
Jx — Jx dynamical model with a small a J [see Fig. 1(a)], 
n[ln(P 2 )] here displays very smooth behavior. Figure 
E^b) depicts a nice linear relation between (In Pa) an d 
ln(iV). The slope of the straight line in Fig. EJb) gives 
the fractal dimension D2 — 0.279, a value close to that in 
the J x -J x model with 0.95 < a J < 1.05. As in Fig.[TJc), 
Fig. [SJc) shows the fitting of the variance of ln(P 2 ) with 
N, using Eq. ([3]). Interestingly, the values of the fitting 
parameters are found to be cr 2 (oo) ~ 0.60, A = 1.33, both 
are similar to those determined in Fig. 1(c). More inter- 
estingly, this fitting shows that [cr 2 (oo) — er 2 (7V)] scales 
as N 1 , with 7 = D^. This supports our finding in 
Fig. QIc) and Fig. 2(c). We have also studied other 
cases of PRBUM-COE using other PRBM as the input 
of Mezzadri's algorithm. For example, we find that if 
the parameter b is set at ~ 0.08, then the D2 of the 



C. PRBUM-CUE 

This ensemble aims to model a critical Floquet opera- 
tor without time-reversal symmetry. All eigenstates of an 
ensemble of PRBUM-CUE matrices are used for our sta- 
tistical analysis. For b = 0.1, Fig. [7{a) displays n[ln(P 2 )] 
versus ln(P 2 ), showing again a smooth dependence. Fig- 
ure Efb) shows the corresponding (ln(P 2 )) versus ln(AT), 
which yields the fractal dimension D 2 — 0.251. In Fig. 
[He), we fit the dependence of ln[cr 2 (cx))-cr 2 (A)] in ln(JV), 
yielding [ct 2 (oo) - a 2 (N)} - N^, with 7 = P 2 /2 [in- 
stead of £> 2 /4 predicted by Eq. (4)]. This also confirms 
our early observations in the J x — J y model. The values 
of the fitting parameters are found to be <t 2 (oo) ~ 0.85 
and A ~ 1.05, which are close to what we found in Fig. 
[3fc). We have also checked that if we perform analo- 
gous calculations for b ~ 0.07, then the D 2 value for the 
PRBUM-CUE ensemble will be close to that found in 
Fig. 4(b). Given these results, we are led to the con- 
clusion that PRBUM as proposed above do share some 
general aspects with periodically driven systems having 
critical eigenstate statistics. 
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FIG. 7: (Color online) (a) Distribution of ln(P 2 ) ob- 
tained for PRBUM-CUE, with the matrix dimension N = 
201(4000), 401(2000), 801(1000), 1601(500) and 3201(250). 
The numbers in the brackets give the size of the ensemble. 

(b) Same as in Fig. GJb) and Fig. 4(b), yielding D 2 ~ 0.251. 

(c) Same as in Fig. EJc) and Fig. 4(c), the fitting gives 
(t 2 (oo) ~ 0.85, A ~ 1.05, and 7 = D 2 /2. 



IV. FLOQUET EIGENSTATE STATISTICS OF 
THE KICKED TOP MODEL 



Finally, as a numerical "control" experiment, we study 
the Floquet eigenstate statistics of the standard kicked 
top model. This will help appreciate the difference be- 
tween a normal driven system and a critical driven sys- 
tem, both of which can have a chaotic classical limit. 
Consider then the following Floquet operator for the 
standard kicked top model |15| . 



Pktm = exp 



■iaJ x ), 



(12) 



which is just the last two factors of Eq. (jfjj), with the 
same parity symmetry and time-reversal symmetry as 
the J x — J x model. In addition, we set the parameter 
r]/J = h v at the same value as given in Eq. ©. We 
construct a statistical ensemble by considering a range 
of a, i.e. 0.95 < a < 1.05 (with chaotic classical lim- 
its). We carry out the Floquet eigenstate statistics in 
the odd-parity subspace, whose dimension is J. Because 
the classical limit is found to be chaotic, we compare the 
statistics with that associated with Dyson's COE matri- 
ces in random matrix theory (RMT). 

Figure^a) and (b) compare II[ln(p2)] associated with 
Pktm with that obtained from COE matrices, for differ- 
ent J. The difference between the actual dynamical sys- 
tem and the COE can hardly be seen. FigureEJc) depicts 
(ln(P2)) as a function of ln(J), with the results of the 
kicked top (open circles) almost on top of those of COE 
matrices (crosses). The solid line in Fig. [5Jc) represents 



the theoretical curve for (ln(P2)) obtained from RMT, 
i.e., (In(p2)) ~ In 3 — ln(J). The agreement between 
numerical COE results, analytical RMT result, and the 
kicked top system as a classically chaotic dynamical sys- 
tem is almost perfect. From the curve shown in Fig. 
IHKc), it is clear that D2 here is unity and as such the 
system does not show critical behavior. This non-critical 
behavior indicates that the Floquet states of the kicked 
top model are essentially random states, a feature fun- 
damentally different from our double-kicked top system 
that has a butterfly spectrum and critical statistics in the 
Floquet cigcnstates. It is also interesting to note that in 
Fig. [8ja) and (b), as J increases, n[ln(P 2 )] becomes nar- 
rower and develops higher peaks. This is an indication 
that, unlike the critical cases studied above, II(lnP2) for 
the standard kicked top model approaches a Dirac-delta 
type singular function with zero width (i.e. <7 2 (oo) — > 0) 
as J increases. 
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FIG. 8: (Color online) (a) Distributions of ln(P 2 ) for 
the standard classically chaotic kicked top model, for 
J = 100(4000), 200(2000), 400(1000), 800(500) and 1600(250). 
The numbers in the brackets are the size of the Floquet matrix 
ensemble. In constructing the ensembles we have considered 
a range of a, i.e., 0.95 < a < 1.05. (b) Distributions of ln(P 2 ) 
for the standard Dyson's COE matrices, with the same matrix 
dimension as in the kicked top model and the same ensemble 
size, (c) Analogous to Fig. [Tfb) and Fig. 7(b), the scaling 
behavior of (ln(P 2 )) vs ln(J) is shown. Open circles are nu- 
merical results for the kicked-top model, crosses are numerical 
results associated with COE random matrices, and the solid 
curve represents the theoretical prediction from the random 
matrix theory. The scaling shows that D 2 = 1 in the stan- 
dard kicked top model, which is dramatically different from 
our observations made from the double-kicked top model. 
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V. CONCLUDING REMARKS 

In this numerical study we have examined the statistics 
of the Floquet eigenstates of a recently proposed double- 
kicked top model with multifractal Floquet spectra. Fol- 
lowing the methodologies used in studies of Anderson 
transition in time-independent systems, we have shown 
that the Floquet eigenstates associated with multifractal 
Floquet spectra also display critical behavior. In partic- 
ular, we focus on the distribution of ln(P 2 ) and examine 
how the quantity (ln(P2)) averaged over all states scales 
with the dimension of the Hilbert space N. It is shown 
that (ln(P2)) scales linearly with ln(iV), with the slope of 
this linear scaling giving the fractal dimension D2 of the 
Floquet eigenstates. The values of D2 are found to be 
far from unity (as a comparison, we showed that similar 
analysis for a standard kicked top with a chaotic clas- 
sical limit yields D2 — 1), constituting strong evidence 
that the Floquet eigenstates are fractal and hence ly- 
ing between localized and delocalized states. Though we 
have worked on P2 only, we note that similar analysis can 
be done for P q defined in Eq. ([1]). One may then define 
a generalized fractal dimension D q and further establish 
the multifractal nature of the Floquet eigenstates. 

The variance of ln(P 2 ), denoted a 2 (N) for a Hilbert 
space of dimension N, is also examined. In Anderson- 
transition studies with PRBM, cr 2 (N) is known to scale 
as N^ 1 with 7 = D 2 /(2/3) for one-dimensional systems, 
where /3 = 1 (2) for a system with (without) time-reversal 
symmetry. By contrast, in our critical driven system, 
<J 2 (N) is seen to scale similarly, but with 7 = D%/ p. This 
reflects an interesting difference between time-dependent 
systems and time-independent systems. Indeed, eigen- 
states of PRBM are to model those of critical Hermitian 
operators, whereas Floquet eigenstates of a critical driven 
system should be understood in terms of critical unitary 
operators. To justify this understanding, we have intro- 
duced a random unitary matrix ensemble called PRBUM, 
with the variance of the matrix elements of the unitary 
matrices following a power-law distribution. We show 
that the eigenstates of PRBUM share many critical sta- 
tistical features with the double-kicked top model. Most 
important, the variance of ln(P 2 ) of PRBUM does scale 
as N~( D2 /P\ which is the same as in the double-kicked 
top model as a critical driven system. We hence antici- 
pate that this scaling property of the variance of ln(P2 ) 
may be general in critical driven systems. These results 
complement the spectral results in Ref. [l4| and should 
motivate further mathematical and theoretical studies in 
critical driven systems. 
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Appendix A: Mezzadri's algorithm 

This is a simple and numerically stable algorithm to 
generate the CUE matrices from an ensemble of com- 
plex random matrices {Zi}, whose elements are Gaussian 
distributed random numbers with mean zero and vari- 
ance unity. In particular, applying the Gram-Schmmidt 
ortho-normalization method to the columns of an arbi- 
trary complex matrix Zi, one can factorize Zi as: 

Zi = QiRi, (Al) 

where Qi is a unitary matrix and Ri is an invertible 
upper-triangular matrix. One can easily prove that the 
above factorization is not unique. Because of this non- 
uniqueness, the random unitar y m atrices {Qi} are not 
distributed with Haar measure [28[, i.e., the {Qi} matri- 
ces are not uniformly distributed over the space of ran- 
dom unitary matrices. Fortunately, this factorization can 
still be made unique by imposing a constraint on the Ri 
matrices. By some group theoretical arguments, it was 
shown [28] that if one finds a factorization such that the 
elements of main diagonal of Ri become real and strictly 
positive, then {Qi} matrices would be distributed with 
Haar measure and hence form CUE. Following these re- 
sults, the major steps of Mezzadri's algorithm are the fol- 
lowing. First, we start with an N x N complex Gaussian 
random matrix Zi. Second, we factorize Z\ by any stan- 
dard QR— decomposition routine such that Z± = QiRi- 
Third, we create a diagonal matrix 

a ( m r NN \ 
A = diag , 1 , 

Vim I |rjvjv|/ 

where {nz} are the diagonal elements of Ri. As a final 
step, we define P- = A~ 1 P i and Q\ = QiA. By construc- 
tion, the diagonal elements of R[ are always real and 
strictly positive, and as such {Q^} would be distributed 
with Haar measure and can be used to form the desired 
CUE. The symmetric COE matrices can be constructed 
from the CUE matrices in a very simple manner. In par- 
ticular, let U be a member of the CUE generated above, 
then it can be shown that V = UU T will be a mem- 
ber of COE. For the generation of PRBUM advocated in 
this work, we propose to replace Zi in the first step by 
a member in the PRBM ensemble that models Ander- 
son transition. Though there is no mathematical theory 
for our procedure, the uniformly distributed eigenphases 
(not shown here) of our PRBUM ensemble thus gener- 
ated suggest its uniform distribution. 
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